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Finite size scaling for the Schrodinger equation is a systematic approach to calculate the quantum 
critical parameters for a given Hamiltonian. This approach has been shown to give very accurate 
results for critical parameters by using a systematic expansion with global basis-type functions. 
Recently, the finite element method was shown to be a powerful numerical method for ab initio 
electronic structure calculations with a variable real-space resolution. In this work, we demonstrate 
how to obtain quantum critical parameters by combining the finite element method (FEM) with 
<n: finite size scaling (FSS) using different ab initio approximations and exact formulations. The 

critical parameters could be atomic nuclear charges, internuclear distances, electron density, 
disorder, lattice structure, and external fields for stability of atomic, molecular systems and 
quantum phase transitions of extended systems. To illustrate the effectiveness of this approach we 
provide detailed calculations of applying FEM to approximate solutions for the two-electron atom 
with varying nuclear charge; these include Hartree-Fock, density functional theory under the local 
density approximation, and an "exact" formulation using FEM. We then use the FSS approach 
to determine its critical nuclear charge for stability; here, the size of the system is related to the 
00- number of elements used in the calculations. Results prove to be in good agreement with previous 

Slater-basis set calculations and demonstrate that it is possible to combine finite size scaling with 
' (— i ', the finite-element method by using ab initio calculations to obtain quantum critical parameters. 

The combined approach provides a promising first-principles approach to describe quantum phase 
I transitions for materials and extended systems. 

a" 

a 

lE; 

> 

CO 

in 

av 

o: 



PACS numbers: 02.70.Dh 64.60. an 05.70.Jk 82.60.-s 



a kais@purdue.edu 



2 



I. INTRODUCTION 

The theory of finite size scaling (FSS) in statistical mechanics can provide us with numerical methods [1-8] capable 
of obtaining accurate results for infinite systems by studying the corresponding small systems. In the past, basis set 
approximation to many-body systems have been combined with FSS to obtain critical parameters, where calculations 
for these systems were done by expanding the wave function in Slater-type basis sets or Gaussian-type basis functions. 
However, a generalization of larger atomic and molecular systems proved itself difficult using these type of functions 
[9, 10]. 

More recently, we have combined FSS with the finite element method (FEM) in order to obtain critical parameters in 
a potential that admits a finite number of bound states and we obtained excellent results by using both finite-clement 
and Slater-type basis functions [11]. For this potential we also established that the finite-difference (FD) method 
can be combined with FSS. The finite difference method, however, is not a variational method since it approximates 
the operators (i.e. the derivative) and therefore can result in errors of either sign [12]. As a result, we observed non 
monotonic behavior of the calculated critical parameters with respect to the number (N) of grid points used [13]. 
Nevertheless, for a large enough number of grid points N one can recover the correct asymptotic behavior in the large 
N limit. For this particular example we used FD to fourth order and in FEM we used C 1 -continuous basis. 

In this paper, we present a method to combine FEM with FSS for electronic-structure calculations near regions of 
criticality. We note that any method where the wavefunction can be expanded in a complete basis (or set of grid 
points) should also work, provided that in the infinite basis limit (N — > oo) it becomes arbitrarily close to the the 
exact wavefunction. The main advantage of finite element methods, over other global basis methods, resides in the 
fact that matrices obtained by these numerical methods are often banded and tend to be sparse since variables are not 
coupled over arbitrarily large distances. This is ideal for describing extended systems since it facilitates parallelization 
[14]. Moreover, due to the polynomial nature of the basis set most integrals involved can be evaluated easily and 
accurately with a better convergence with respect to the basis set size. 

We examine the two electron system with different levels of approximations to the ground state energy, by solving 
the system using FEM both exactly and using mean-field equations and demonstrate that this formalism is general and 
that it has potential for more complex systems by solving Hartree-Fock equations or Kohn-Sham equations in density- 
functional theory. The paper is organized as follows: In Sec. (II) we discuss FEM formalism, in Sec. (Ill) we discuss 
the FSS approach to quantum systems, in Sec. (Ill) and (IV) we combine FSS with FEM to compute calculations of 
critical parameters for two electron atoms, ,and, lastly in Sec. (V) we discuss the results and applications to extended 
systems. 



II. FINITE ELEMENT METHOD IN ELECTRONIC STRUCTURE CALCULATIONS 

Computational methods in quantum chemistry typically make use of Slater or Gaussian basis sets to approximate 
solutions in atomic or molecular systems. The advantages of these methods are that some of the integrals involved 
in the calculations are analytical and very accurate results can be obtained for these integrals. However, there are 
problems in extrapolating to the complete basis limit. Plane wave methods, having a complete basis, present some 
difficulties if non-periodic boundary conditions are to be used such as in molecule calculations. A common trait of 
all these methods is that they employ global basis, which can result in some difficulties for the formulation of large 
problems. First, the matrix formulation becomes dense using global basis and therefore ill-conditioned for parallelism 
[14]. While for some complex objects, convergence issues may arise if there are large variations of potential terms in 
small regions which cannot be captured by the traditional Gaussian basis set that are centered on the nuclei. 

The FEM presents an alternative solution to these problems. This method uses a variational method to solve 
differential equations by discretizing a continuous solution to a set of values in sub-domains, called elements, and 
employs a local basis expansion in terms of polynomials with variable real-space resolution, allowing for convergence 
to be controlled systematically. The locality of the bases in real space results in sparse and banded matrices, where 
the number of operations for matrix- vector multiplication can be reduced from 0{N Log N) for dense matrices to 
0(N) for sparse matrices, where N is the dimension of the matrix [14]. Some of the first applications of FEM 
in engineering date back to 1950, in problems such as elasticity and structural analysis in civil engineering and 
aeronautical engineering, respectively. In the 1970s, some applications to quantum mechanical problems appeared, 
but were limited to small systems due to the storage limitation [15, 16] Recently, advances in high-performance 
computing make FEM a viable alternative to the traditional approaches of electronic-structure calculations [17, 18]. 
The sparsity of the global matrices resulting from the FEM makes the parallel formulation more affordable. Alizadean 
et al, have recently applied a divide- and- conquer method to the FEM-Hartree-Fock approach for electronic calculations 
showing a facilitation on parallelization and reduced scaling for larger systems [19]. 
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We briefly describe the FEM formulation for a given Hamiltonian to be used in later sections. Let us consider the 
following one-particle Hamiltonian: 

-Iv 2 *i + y* j =e i * i (1) 

where and are the energy eigenvalue and eigenvector respectively, and V is an arbitrary potential that will 
represent a sum of terms such as the electron- nuclear potential (V en ), Hartree potential (Vh) and exchange-correlation 
potentials (V xc ). We can reformulate this problem into a weak form by taking an inner product of an arbitrary "test 
function" $ with Eq (1): 

f $[_I v 2 + v-e 4 ] v i( m = o (2) 

Jn 1 

This relaxes the problem into a variational form where instead of finding an exact solution in the domain f2, we find a 
solution that satisfies Eq.(l) in an average sense everywhere. We then arrive at an equivalent equation by integrating 
by parts (in real space): 

f ^V$(r) • V*i(r) dSl- f J$(r)V*i(r) • h aT + [ $(r)(V - e t ) ^ t (r)dn = (3) 

where T is a boundary term and n is the outward unit normal at r. In the absence of external fluxes we can employ 
Neumann boundary conditions [V^j(r) • n = f(r) r £ T]. Furthermore, in radial coordinates (dfl = r 2 dr) we can 
write the surface term in the following form: 

jf ^*(r)V*i(r) • h dT = [l$(r) f(r) r 2 ]^ (4) 

This surface term vanishes trivially at r = 0, while for r = r c , it can be removed if we choose a test function 
$(r) that vanishes at r = r cut . Since Eq (3) is satisfied for any we can choose $ = <£* and in this manner it 

reduces the problem to an eigenvalue problem. The wavefunction ^i(r) can be then discretized using the Garlekin 
approximation [15, 16]: 

|*(r))«N>)=]T^(rM> (5) 

i 

where Ni(r) are piecewise-linear C° continuous functions (i.e. continuous but not necessarily smooth) and |V>i) arc 
the nodal degrees of freedom. Smoother bases can be constructed by using higher order polynomial (C™ continuous, 
n > 1) and by extending the nodal space to include the derivative of the wavefunction at the node, viz, 1^) = 
J2 i Ni(r)\ipi) + d x Ni(r)\d x ipi}] this ensues continuity of the wavefunction between the nodes. In the nodal basis, the 
generalized eigenvalue problem reads 



J^mHi^j) = e^ilSi^) (6) 

where the above terms can be written in terms of local basis representation, viz 

Hij = [ [^ViV, • VNj + VNiNj]dCl 
JQ 2 

S i:i = [ N t NjdQ. 
Jn 

The finite element method takes advantage of the strict locality of its basis to yield matrices that are sparse 
and structured and due to the polynomial nature of the basis most integrals can be easily evaluated, as in global 
basis approaches. Moreover, since the method is variational, all errors are of the same sign, leading to monotonic 
convergence often from below [20]. 
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III. FINITE SIZE SCALING 



We now discuss the role of finite size scaling in systems exhibiting critical behavior. Phase transitions are associated 
with singularities of the free energy that occur only in the thermodynamic limit [21, 22]. Near criticality, fluctuations 
in the system become correlated over distances of the order of the correlation length £, while for second order phase 
transitions, this correlation length £ diverges. 

The theory of finite size scaling can be used to extract universal parameters describing the transition by studying 
systems of finite size [1, 2, 4-8]. Any physical quantity (P) that would diverge in the thermodynamic limit, becomes 
now "rounded off" in the finite systems (of size L) by a scaling function f(x) [1, 2] 

where t = (T — T c )/T c is a reduced critical parameter describing the phase transition and ^ is the correlation length 
of the system in the thermodynamic limit. 

While classical phase transitions are governed by the fluctuation of temperature that tend to minimize the free 
energy, at zero temperature, these fluctuations are driven by Heisenberg's uncertainty principle [23] . A signature of a 
phase transition is non-analyticity of the ground state energy E(X) with respect to a parameter A in the Hamiltonian. 
The non-analyticity can be the result of avoided crossing, or actual level crossing. In the first case, the configuration 
of the system that minimizes the energy maintains the energy gap (AE\) between the first excited energy state and 
the ground state as the parameter A is varied. In the latter case, the energy levels cross so that the energy is now 
minimized by a configuration that previously corresponded to the first excited state of the system. 

To pin down the transition point involves computing the energy accurately near the critical point. However, this can 
be quite a challenge since most real systems lack an analytical solution to the ground state energy and approximation 
methods, such as perturbation theory or mean-field methods (Hartree-Fock , DFT) which are viable options in stable 
systems, in general are not well defined or become ill-conditioned near critical points. 

This is where the theory of finite size scaling can be recast in the context of the Hilbert space. Any given ap- 
proximation to the ground state energy will result in a "rounding" effect from the "true" critical behavior, since any 
approximation to a wavefunction can be expanded in some basis and truncated to some order N; that is, 



N 

|*> = (8) 

i 

FSS has been used in quantum systems [24-33]. In this approach, the finite size corresponds not to the spatial 
dimension but to the number of elements (N) in a complete basis set used to expand the exact eigenfunction of a 
given Hamiltonian [34-38]. In the variational approach singularities in different mean values will occur only in the 
limit of infinite number basis function, viz., 

<0>r~(A-A c )^ (9) 

The FSS ansatz in the Hilbert space takes the following form (for expectation values of operators acting on the 
Hilbert space). There is a scaling function F(x) such that 

(0){ N) ~{0)?F (N\\-\ c n (10) 

where A is a critical parameter in the Hamiltonian, O is a quantum mechanical operator near criticality, and Fo is a 
scaling function for each different operator but with a unique scaling exponent. At the critical point, the expectation 
value is related to TV as a power law, (0) (A ° - N-^'". In this manner the scaling function Fo will remove the 
divergence as it is computed in the finite basis. 

The asymptotic behavior of operators near the critical point, in this case the ground state energy, can be expressed 
using critical exponents similarly to classical phase transitions, viz 



A£ A ~(A-A c r (11) 

where a is a universal critical exponent. In addition to the vanishing energy scale, a diverging length scale also occurs 
which is associated with the exponential decay of equal-time correlations in the ground state of the system [23]. For 
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our case, we focus on the behavior of the wavefunction at threshold, where the correlation length is defined by: 

e ~ |A - A c |— (12) 

In order to apply finite size scaling method to quantum system, we assume a Hamiltonian of the form 

H = H + V X (13) 

where A is the coupling parameter showing a phase transition at A = A c . For potentials linear in A ( V\ = XV), 
Simon [39] showed that the critical exponent is equal to one if and only if the Hamiltonian H(\ c ) has a normalizablc 
eigenfunction with eigenvalue equal to zero. The existence or absence of a bound state at the critical point is related 
to the type of singularity in the energy. Using statistical mechanics terminology, we can associate a "first order phase 
transition" with the existence of a normalizable eigenfunction at the critical point. The absence of such a function 
could be related to a "continuous phase transition." To obtain the critical parameters, we define the following function 
that will remove the N-dependence at a critical point: 



A (X;N,N') = 



H(Q)x/(Q)x') 



(14) 



]n(N'/N) 

For the energy operator O = H, and using the customary a Greek letter for the corresponding exponent [i Q we have 



A H (A C ;A,A') = - 



(15) 



To arrive at an expression that isolates the critical exponent a from numerical calculations, we make use of the 
Hellmann-Feynman theorem [40] 



dX \ dX 

At this point it becomes convenient to define another function that will be used in the following sections. 



T a (X,N,N') = 



A H (X;N, N') 



A H (X;N,N')-Ae^(X;N,N') 



(16) 



(17) 



At the critical point this function is independent of N and N' and takes the value of a. Namely, for A = A c and any 
values of N and N' we have 



T a (X c ,N,N') = a 



(18) 



IV. IMPLEMENTATION 



In the following, we illustrate the combined approach (FEM+FSS) for the two-electron atom, but emphasize that 
the method is general and can be applied to more complex atomic and molecular systems. We seek the minimum 
charge (critical charge) needed to bind both electrons to the nucleus with charge Z. It is well known that ignoring 
the electron-electron energy correlation will not predict the stability of H~ [41], since the critical charge obtained 
by such assumption is Z c w 1.02. Stillinger (1974) [42] took into account the e-e correlation by using a non-linear 
wavefunction e -<wi-&r 2 +e -an-6r 2 in order t0 find the 

minimum energy as the nuclear charge Z was varied. He found 
a critical charge around Z c rts .9537. Baker et al (1990) [43] performed a thorough 401-order perturbation calculation 
to resolve the controversy for the critical charge; they obtained a critical charge of Z c w .91103. This is the value 
we use as reference. In the following sections wc illustrate the mean-field equations and exact formulations for the 
ground state of the two-electron system, while in the next section we apply FSS to these approximations. Let us start 
with the Hamiltonian for the two-electron atoms; it is given by 

H = - \vl + - (1 + -), (19) 

2 2 \r 1 -r 2 \ n r 2 
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where r\ and r 2 are the coordinates of the two electrons referenced from the nuclei. We use Hartree units (h = m = 
e = 1) throughout in our calculations. 

Hartree Fock (HF) Approximation: 

In the Hartree approximation the total wavefunction is approximated as a product of single particle wavefunctions 
is ^{r\,r 2 ) — VrfVi)^^)- We can arrive at a self-consistent form by multiplying the Hamiltonian by the complex 
conjugate of the aforementioned wavefunction and integrating one coordinate out, we arrive at 



- ^ + / d*r 2 -±— |V>2(r 2 )| 2 ]Vi(ri) = 6^(n) 

z ' 1 J | ' 1 '2| 



(20) 



Using the FEM (see Sec. II) this expression can be reduced to a generalized eigenvalue problem H k iipi = eS^i- The 
Hamiltonian H k i is written in terms of the single electron energy term h k i (i.e. — ^V 2 — ^) and a direct coulomb 
term J k i defined below, viz 

h J} Jm 

ri — : c— 

H kl 



= J [\ VN k * VTVf + ^N k Nf + g Qmj^A^ri (21) 



Qklij^NkinYNiin) I ^_ r ^ N t (r 2 yN 3 (r 2 )d A r 2 



Sm = J N k N?d 3 n 

Since electrons are fermions, the product of the single particle states needs to be antisymmetric with respect to 
exchange of the coordinates on the electrons. This can be achieved by writing the wave function as a Slater determinant 

[44]: 



*(r-i,r 2 ) = 



^i(r 2 ) 4> 2 (r 2 ) 



(22) 



Searching for the best Slater determinant self-consistently in the Hamiltonian [Eq.19] leads to the Hatree-Fock 
equations [44-46]. 



/ 



r i (ri)(-\vl---e)Mri)d 3 r 1 + \ f f -J— |^ 2 (r 2 )| 2 -^(r 2 )V 2 *(r 2 )^^Vi(r-i)V2(ri)}d 3 r 1( i 3 r 2 = 
2 n 2 J J n - r 2 n - r 2 \ 



(23) 

where the additional term on the right is known as the exchange Coulomb term, denoted by K k i , and tends to repel 
states of the same (parallel) spin. We can write this term in the local basis form by looking at the direct Coulomb 
term and noting that the indices € ip 2 and similarly indices (k, I) G ipi. The operator K M is obtained from J ki 
by switching indices (j -H> I), which effectively switches the coordinates n and r 2 on two of the wavefunctions and 
makes this operator non-local. 

The ground state energy (for N e electrons) is determined by solving for the ground state eigenfunction \tp g ) for the 
eigenvalue problem (Eq.23) and computing the expectation value below (to avoid double-counting terms), 

Ehf = J> fl | hj - - (Jj -K 3 )\iP g ) (24) 

3 

Note that for our case, the exchange term does not arise since the ground state of the helium atom does not have 
parallel spins. In order to improve upon the Hartree-Fock approximation wc can take into account the correlation 
energy from an approximation to the free electron gas due to Wigner. Here the "correlation" potential can be 
obtained from Refs. [47, 48]. Note that the references cited use Rydbergs units, and reference [47] additionally uses 
e 2 = 2, hi 2m = 1. In Hartrees the correlation potential is given by: 
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"■< r > - < 25 » 



where r s is the local density parameter for the atom, that is 

P(r) = ~4~~~~3 (26) 

P(r) = P(rh + p{r)i (27) 

The last expression follows for helium since in the ground state the two electrons have opposite spin, but share the 
same spatial orbital. In this manner the correlation potential is included self-consistently and the density is obtained 
from the sum of the densities for each electron i.e., p(r) = 2\tp(r)\ 2 . Finally, the correlation energy is obtained as 

E C orr = (lpg\Vc(r)\lp g ) 

and the value for the total energy is given by: 

Etotai — E H f + E corr (28) 



Density Functional Theory 



In the Hartree-Fock approximation, the many-body wavefunction is written in the form of a Slater determinant 
while the correlation energy is treated separately. This can be done for small systems, but for larger systems this is 
not an easy task. The theory of Kohn and Sham developed in 1965 simplifies both the exchange and correlation effect 
of electrons by reducing their effect to a effective potential Vks that is in principle exact. In their seminal paper [49], 
they showed that to determine the ground state energy it is sufficient to determine the ground state electron density 
p(r). By constructing a simple system of non-interacting electrons which yields the same density as the "real" fully 
interacting system, the ground state energy can be obtained. 

The Kohn-Sham equations for two electrons are: 



~V? - ^ + J <* 3 r 2 * P(r 2 ) + V xc (p)]Vi(ri) - «M^i) 



(29) 



These equations are self-consistent in the density and differ from the Hartree-Fock equations only in the exchange- 
correlation potential V xc [p] which is a function of the total local density: p(r) — J^f I^M 7 ")! 2 

V xc [ P {r)\ = j^E xc [p(r)} (30) 

While the exact form of the exchange-correlation functional E xc [p] is not known, the simplest approximation is the 
local density approximation (LDA) 

E^ A [p(r)]=E x [p(r)]+E c [p(r)} (31) 
E x [p(r)]=-^)(^Jp(r)^d\ (32) 



It follows from Eq.(30) that the potentials are given by: 



V x [p(r)} = -(^) 1 /3 /9(r) V3 (34) 

7T 
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TABLE I. Comparison of various components to the Kohn-Sham and Hartree-Fock energies in Hartrees for helium. The 
reference values (Perdue-Wang, Exact) are correct to all digits quoted [50, 51], while the accuracy of the FEM estimates 
on Etot is expected to be on the order of two digits (0.01) for N=100 and three digits (0.003) for N=200 for both LDA and 
Total Energy, refer to the discussion on the text for details. 







LDA 




Total Enerj 






FEM 


FEM 


Perdue- Wang 


FEM 


FEM 


Exact 




(N=100) 


(N=200) 


[50, 51] 


(N=100) 


(N=200) 


[50] 


Etot — Ekitl + -E'en + Er + E xc 


-2.813 661 


-2.821 852 


-2.834 455 


-2.886 297 


-2.904 925 


-2.903 724 


Exin 


2.703 352 


2.726 064 


2.767 389 


2.806 848 


2.855 856 


2.867 082 




-6.538 681 


-6.571 185 


-6.624 884 


-6.664 319 


-6.737 783 


-6.753 267 


Eh 


1.973 543 


1.978 456 


1.995 861 


1.020 329 


1.026 754 


1.024 568 


E x 


-0.850 814 


-0.854 086 


-0.861 740 


0.000 000 


0.000 000 


0.000 000 


E c 


-0.101 061 


-0.101 092 


-0.111 080 


-0.049 135 


-0.049 814 


-0.042 107 


e 


-0.563 621 


-0.565 805 


-0.570 256 


-0.957 555 


-0.963 991 


-0.903 724 



T/r^l ~°' 29 1 0-29r s (r) 



and the ground state energy of the system is obtained from [44] : 



E DFT =it ei ~\j J rfrtcPnpin) ]_ p(r 2 ) + E xc [p] - J d 3 rV xc [p(r)] (36) 

i 

The Kohn-Sham equations are in general simpler than HF equations, while including correlation effects beyond the 
HF approximation. Furthermore their computing time grows as a power of the number of the N-electrons, whereas 
the exact solution for the N-electron Schrodinger equations is exponential in N [46] . 



Let us comment on the difference between the two approaches, even though both Hartree-Fock and Kohn-Sham 
equations have similar potential terms: Vkin ,V en ,VH ,V XC [see Eqs.(23) and (29)]. In HF the energy is calculated 
using the wavefunction explicitly, whereas in DFT the energy is a functional of the density. For example, for the 
correlation energy we have in the HF approximation: 

E* F = J \yj(r)\ 2 V c (r)dr 3 

while for the DFT approximation we have: 

E? FT = J p(r)V c [p(r)]dr 3 

Care should be exercised when computing the LDA exchange potential since this term is non-linear in the density 
and as a result the normalization does not trivially cancel out when solving the eigenvalue problem. For example, 
in this work we use the normalization J \ijj(r)\ 2 r 2 dr = 1 so that a factor of y^(47r) is absorbed into the wavefunction 
ip(r) — > ip(r)/ y/lAw). This in turn changes the density to p = 2\ip(r)\ 2 — > 2ip(r) 2 /4-k, and therefore the exchange- 
potential in LDA becomes 



In order to satisfy the boundary conditions ip(r) — > as r — > oo, we use a uniform mesh with a cutoff radius 
at r cut = 10 a.u. The radial equations where solved using C°-continuous polynomials, that is, Nl^,r(x) = x and 
Nl^r(x) = 1 — x; where the function Nl^.r(x) interpolates the wavefunction from the left node to the right node 
and similarly for Nl^r(x). We solved all the integrals numerically using 10-point Gaussian quadratures and in order 
to solve the generalized eigenvalue problem we have used a fortran subroutine EWEVGE [52]. 

In Table I we show the energy components of the Kohn-Sham equations under LDA approximation and the Hartree- 



FIG. 1. example of an adaptive mesh [56] used in approximating the solution in three dimensions. In this case the number 
of elements is N e i ement s = 15 x 15 x 3 = 675, while the size of the element progressively increases by a factor of 1.3 as ri,2 is 
increased 



Fock equations using the Wigner correlation potential for the helium atom. We note that the reference values for 
the Total Energy (HF+Wigner correlation) are exact, therefore even in the infinite basis limit the expression will not 
converge to the value shown. We also point out that the Perdue- Wang values on the correlation energy are obtained 
by using the random-phase approximation (RPA) for the free electron gas [51]. Our convergence value is be limited 
by how well we can approximate this value. We have checked our level of convergence by increasing the number 
of elements to N=1000 in a linear mesh of size r cut = 10 a.u with C°-continuous basis (linear polynomials). We 
obtain a total energy of Ef Q fJf = 2.824596 and correlation energy of E^ EM = —0.101103, and when compared to 
the values by Perdue- Wang (see Table I ), the errors are 5E to tai — 0.009859 and SE C = 0.009977 respectively. We see 
that the error in the total energy is dominated by the correlation term arising from using Wigner correlation rather 
than RPA description. However, taking this difference into account we can get an estimate for the error of all other 
terms, viz 8E Kin + SE en + 5E H + SE X w SE to tal - &E C = 0.000118. In FEM, the error on L 2 integrals such as the 
energy is 0(h^ p+1 ^), where h is the element width and p is the order of the polynomial used. Thus, for our example 
using N=1000 in a uniform grid size with r cut — 10 a.u. we have h = 0.01 , and since we employ a linear basis, the 
error we expect is 0(h^ 1+1 ^) = 0.0001. In this manner, if we take into account the intrinsic error due to the different 
energy-correlation terms we can recover the level of convergence expected by using FEM. 



Exact Formulation 



We can compare the one-particle formalism above with an exact formulation of the system. In 1930 Breit showed 
that the Schrodingcr equation for the two-electron atom becomes separable if one of the coordinate axes is aligned 
with one of the radius vector [53]. For the s state, the exact wavefunction depending on six total variables *(r"i,r 2 ) 
can be reduced to only three variables: the positions of the two electrons r\ and r 2 , and the angle between them 12 , 
viz 

10**0* 10**0* 1,1 1.0* 0* T „, Z Z 1 _._ 2 2j , , n . . 

o^— ^- + o^— ^-+o ("2+"2 )« TTa ^ + * + £)*r?r2dndr 2 dcos0i2 = 37 

2 0n 0ri 2 or 2 or 2 2 rf rs #cos#i 2 0cos#i 2 n r\ r i2 



where 



r i2 = \Jr\+rl- r x r 2 cos6>i 2 



We approximate the solution to the above Hamiltonian using FEM by discretizing each variable above independently 
as in [54]. The boundary conditions implemented make the wavefunction vanish at a cutoff radius of ri ;2 = 40 a.u, 
and the angular range is between — 1 < cos #12 < 1. Furthermore an adaptive grid was refined closer to the nuclei 
as shown in Fig.(l) to speed up calculations. In our implementation we use C 1 -continuous basis, that is Hermitc 
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polynomials Nl^r(x) = 1 — 3x 2 + 2x 3 and Nl^^x) = 3x 2 — 2x 3 to interpolate the left and right nodal values 
of the wavefunction and additionally Nl^h(x) = h x (x — 2x 2 + x 3 ) and Nl^r,( x ) = h x (x 3 — x 2 ) to interpolate the 
derivatives of the wavefunction, where h x is the element width and x — {ri, r 2 , cos9i 2 }. In this manner, the continuity 
of the wavefunction is enforced by using this basis. All the integrals involved in constructing the matrices were solved 
numerically using 10-point Gaussian quadratures, and in order to solve the generalized eigenvalue we used a sparse 
matrix solver in ARPACK++ [55]. 

For the mesh shown in Fig. 1, we obtain a ground state energy of E® xact — —2.7578 using C°-continuous basis and 
E\ xact = —2.8994 using C 1 -continuous basis. The latter value is close to the energy obtained by Levin using the same 
approach = -2.9032 [54], which is very close given that we perform all integrals numerically. 



V. RESULTS AND DISCUSSION 



In order to find the critical parameter we approximate the energy gap in the ground state numerically as 

AE (Z)^E 2e (Z)-E { e \z), 

where E^ e ^ is the binding energy obtained for two electrons and E^ is binding energy for a single electron. The 
theory of finite size scaling is applied to a helium-like Hamiltonian of the form, 

H = -hi \v\ + -(- + -), (38) 
2 2 \n-r 2 \ n r 2 
. ' v . ' 

H V z 

where r\ and r 2 are the coordinates of the two electrons referenced from the nuclei, and Z is the charge of the nucleus, 
this will play the role of the critical parameter. For reference, we note that the transformation r — > r/Z takes Eq. 
(38) into: 



H(r - r/Z) = Z»[- * V? - \v 2 - (I + I) + 1^] (39) 
v . ' * . ' 

H' a V' z 

this form is often used in the literature to obtain the critical parameter A c = 1/Z C (see [42, 43]). Below we outline 
a general procedure that can be used to combine FSS with numerical methods that expand a given wavefunction in 
the Hilbert space in any complete basis at a truncated level N. 



General Procedure: 



1. Compute the expectation value of (Hy ' and (^if") f° r various values of N (number of elements). In this 
case the Hellmann-Feynman theorem (see Eq. 38) becomes useful since in this case Hz is linear in Z and the 

derivative is quite simple {^f-) — >• (^f") 

2. Construct the gamma function T a (Z, N, N') (see Eq. 17) for several values of (N,N') vs Z. At criticality, the 
gamma function T a (Z,N,N') becomes independent of N and N' and the curves will cross. In Fig. 2 we show 
the gamma function we obtain by using different number of elements and how it begins to cross for various 
approximations to the energy. 

3. In order to go to the complete basis set limit {H)°° , we have to let (N,N') — > oo. The best choice is to keep 
N — N' minimal [1, 6]. Then, we proceed as follows: We solve (Hy' for (N-l, N , N+l). This provides us 
with two gamma curves: T a (Z 7 N — 1, N) and T a (Z, N, N + 1), with minimum N — N' . The Crossing of these 
two curves gives us the pseudo-critical parameters and Z^ with a truncated basis set number N. Fig. 3 
illustrates the N dependence of the crossing. 

4. In order to extrapolate the values and Z^ to the infinite N limit, we use the algorithm of Burlish and 
Stoer [57]. Fig.s 4 and 5 show their behavior by systematically increasing N for the mean- field approximations 
and exact formulation, respectively. 



11 




1.028 1.03 1.032 1.034 1.036 1.038 1.04 1.042 1.044 



z 

(a) 




0.91 0.911 0.912 0.913 0.914 

z 



1.015 
1.01 



r(145,150) 
T(1 95,200) 
T(295,300) 
T(395,400) 
T(495,500) 
r(595,600) 
T(695,700) 



8 1 



005 



0.995 



0.99 




0.91 



0.92 



1.2 



1.15 



1.1 



1.05 



0.95 



0.93 



0.94 



0.95 



0.96 



z 

(c) 



T(1 92,243) 
r(243,300) 
r(300,363) 
T(363,432) 
T(432,507) 
r(507,588) 
r(588,675) 
T(675,768) 
T(768,867) 




0.915 



0.916 



0.917 0.91 i 

z 

(d) 



0.919 



0.92 



FIG. 2. (Color online) Plot of r a (Eq. 17) obtained by using FSS method as a function of the charge Z for increasing number 
of elements (N) for: (a) Hartree-Fock only, (b) Total Energy (HF + correlation energy), (c) LDA approximation, and (d) exact 
formulation 



While finding values and Z^ N ^ and extrapolating to the l/N — > limit some numerical fluctuations can 

arise depending on the approximation used. We discuss some ways of dealing with these fluctuations. First, one 
can improve the basis used from C°-continuous polynomials to (^-continuous polynomials in order to obtain greater 
accuracy in the calculations, however, matrices double in size increasing the memory cost as well as the computing 
time. Second, one can construct a gamma function from a potential in the Hamiltonian that explicitly takes into 
account all variables. For example in the three dimensional (3D) case we found that it is better to use the transformed 
Hamiltonian (Eq.39) where the perturbation V\ is given by V\iz — i , , rather than Vz = — + — • This 

does not change the physics of the system since both expressions can be used to find the critical value: either 1 jZ c or 
Z c , which corresponds to one of the electrons going from a bound state into the continuum. The first perturbation 
V\/z takes into account all three variables when the expectation value is computed whereas the latter Vz introduces 
numerical errors in the expectation value. 

Last, another method is to increase the minimum difference \N — N'\ > 5 used in finding a crossing between 
T a (Z, N — S, N) and T a (Z, N,N + 5) until a systematic crossing is obtained for all values of N; in the LDA approxi- 
mations we used \N — N'\ = 5, since the Gamma curves change very little as the number of elements increases (see 
Fig.2-c), whereas in the 3D formulation we used \y/N/3- y^VV3| = 1 since in this case the number of elements 



12 



1.03 
1.025 

1.02 
1.015 



1.005 



0.995 



0.99 



r(99,100) 

r(ioo,ior 

r(1 99,200) 

r(200,201 

r(499,500) 

r(500,5or 



1 - 




0.91 



0.915 



0.92 



0.925 



0.93 



FIG. 3. The crossing of T a (Z, N — 1, TV) and T a (Z, N, N + 1) as N is increased using Total Energy approximation in HF. The 
crossing has been slightly enhanced in order to demonstrate the dependence of N. 



TABLE II. Critical parameters extrapolated from Fig. (4) and Fig. (5) 





HF 


LDA 


Total Energy 


Exact 


Reference [10, 43] 


critical charge 












z c 


1.03114 


0.92808 


0.90946 


0.91857 


0.91103 


critical exponent 












a 


0.99971 


1.00011 


1.00009 


1.00082 


1.00000 



increases as N e i ements — N ri x N r2 x N cos g, and N cos g was kept fixed at three. In solving the self-consistent equations 
[ Eqs.(23) and (29)], we found that the region below Z < 1 became problematic when updating the self-consistent 
equations. This is because as the nuclear charge decreases, the wavefunction decay becomes slower with respect to 
r and the cutoff radius might lead to unphysical solutions (i.e. solutions which depend on r cut ). To get around this 
problem, we found it is useful to perturb the wavefunction with the arithmetic average of the current and predicted 
wavefunctions, rather than just the predicted wavefunction. In this manner, we can avoid converging to a local minima 
that corresponds to an unphysical solution. 

In phase transition theory it is well known that using mean field theory (for dimensionality < 4) does not give 
correct values near a critical point. This is indeed the case by using only the Hartree-Fock approximation (see Table 
II ), however we see that incorporating the exchange-correlation energy within the LDA approximation improves 
substantially the critical charge within 2% of the accepted value and within 1% using Total Energy approximation 
and exact formulation. It is interesting to note that all approaches give a consistent estimate for the critical exponent, 
regardless of the approximation used. This study demonstrates that Finite-element method is a viable tool that can 
be applied in combination with finite size scaling in order to obtain good estimates of critical parameters starting 
with ab initio calculations for electronic calculations such as DFT and HF and compare well with other methods that 
are deemed more precise. 

We now check the consistency of our approach. In the theory of continuous phase transitions, a quantity diverging 
in the thermodynamic limit (N — > oo), becomes a regular function for finite N, viz (0) {N) ~ N-< i "' v . We can check 
that this assumption is indeed satisfied for any value N. In the infinite basis limit we have the following form near 
criticality: (O) (oo) - (Z - Z c y° (Eq.ll). In this manner, can obtain the m th derivative of the energy near criticality 
given that we have a limiting expression for the energy, viz 



E z ~ (Z - Z c ) a 



d m E z 
dZ m 



(Z - Z c ) a ~ 



(40) 
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FIG. 4. Extrapolated values for the critical exponent a (top) and Z c (bottom) for: (a) HF approximation, (b) LDA approxi- 
mation, and (c) Total energy. The filled dots are the extrapolated values for the 1/N — > limit 
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FIG. 5. Extrapolated values for the critical exponent Z c and a in the 1/N — > limit solving the system exactly in 3D. We 
construct the F a function using Hermite interpolation polynomials for two potentials: Vz = —Z(^- + (denoted by □) and 

(denoted by Q)j the latter is obtained by making the transformation r — > r/Z in the Hamiltonian. 



V z = 



XIZ 



\-r^ — r\T2 cos 6 

We checked that the fluctuations seen in the figure vanish by replacing the numerical approximation to the binding energy with 
the analytical form, — Z 2 /2, while the extrapolated value remains unchanged. 



In a finite basis the above expression becomes a power law in the number of basis with the following characteristic 
exponents. 



E { Z N) - N~ a ' w and 



d m E 



(N) 



dZ r - 



(41) 
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Excluding the trivial case N=0, one can therefore use Taylor expansion around \Z — Z c \ < 1. 

E iN) T <r^P ( Z -z c r m 



= V N ( - a+m)/v { - c) (43) 

m! v ; 

m 

= jV -^ E [^^-^r (44) 

^ m! 



= N- a ' v G (N 1 / v (Z - Z c )) (45) 

In this manner we obtain a universal scaling function Go- Plotting E^ N ^N a / v vs the argument of the function 
Go would confirm whether the scaling function depends on N or not. Fig. 6 shows that all quantities obtained from 
different N values do in fact collapse into one curve, and this confirms the scaling hypothesis for FSS with FEM 
(Eq.10). In order to make the data collapse, it should be emphasized that the values a and Z c used to produce 
Fig. 6 correspond to those obtained in Table II. If other values are used, the curves will not collapse and they will 
be distributed in a family of curves (N dependent) all deviating away from the unique curve shown in Fig. 6. The 
correlation length exponent was varied until the best collapse curve was found. This occurred around v rts .85 ± .05 
for all implementations. 

The solution to a set of Hartree-Fock or Kohn-Sham equations for large number of atoms becomes the limiting 
step, and it is here that the FEM can extend the range of sizes that can be investigated by other current means. 
The efficiency of FEM in the context of large-scale calculations derives from the strict locality of the FEM basis in 
real space, thus minimizing the need for extensive interprocessor communications on massively parallel computational 
architectures [14, 20]. The work of Pask et al [14, 20] has formulated a real-space DFT approach for condensed 
matter applications without the need of using expensive Fourier transforms. Impressive results are obtained for 
GaAs calculations using this approach on DFT+FEM that systematically converges, as the number of elements is 
increased, to the values obtained using the conventional plane wave calculations, employing Fourier transforms and 
Ewald summations. 

The scope of the dimension attainable with FEM has also been demonstrated by Sterne et al [58]. They have 
employed FEM to compute positron distributions and lifetimes for various defect structures. Because the distribution 
is non-electron- like (concentrating away from the nuclei rather than around them), the generality of the FEM basis 
makes it particularly well suited for such problems. In order to study defects they have modeled a 5488 atom cell; as 
far as we know, this is one of the largest such calculations reported by any method to date. 

Another applications of FEM to large-scale ab initio calculations is the work of Tsuchida et al on molecular 
dynamics simulations of liquids [59]. They apply FEM study of liquid formamide (HCONH2) using a generalized 
gradient approximation to the exchange and correlation terms to improve the modeling of hydrogen bonds. The use 
of adaptive coordinate transformations helped to double the resolution of the basis at the C, N and O centers. The 
liquid was modeled by 100 molecules in a cubic superccll of side 35.5 a.u., for a total of 600 atoms. The scope of 
this simulation is of the order of the largest one performed by standard PW methods to date and demonstrates the 
capacity of FEM for such large-scale computations already. 

A new promising approach developed recently by Mazziotti [60] is the large-scale semidefinite programming for 
many-electron quantum mechanics. In this approach, the energy of a many-electron quantum system can be ap- 
proximated by a constrained optimization of the two-electron reduced density matrix (2-RDM) that is solvable in 
polynomial time by semidefinite programming (SDP). The developed SDP method for computing strongly corre- 
lated 2-RDMs is 10-20 times faster than previous methods [61]. This approach was illustrated for metal-to-insulator 
transition of H 50 , with about 10 7 variables. 

These examples along with our current work, suggest that we can combine the Finite-element method with finite 
size scaling to obtain critical parameters in extended systems. One possibility is to calculate the phase transition at 
zero temperature in the 2-dimensional electron gas. The ground state energy of the system favors a crystallize state 
as the density of the electrons is lowered. Tanatar and Ccpcrey [62] performed a variational Monte Carlo simulations 
to the total energy using Hartree-Fock approach 

E (r s ) = E HF (r s ) + E c (r s ) (46) 



They find the gas-to-crystal transition at r s = 37 ± 5. This system would be ideal to test FEM with finite size scaling 
for quantum criticality in extended systems. 
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(c) (d) 

FIG. 6. (Color online) Data collapse for various values of N using: (a) HF, (b) Total Energy, (c) LDA approximation, and (d) 
exact formulation. The value from this collapse study that give us the best correlation length exponent v ~ .85 ± .05, it should 
be emphasized that the values a and Z c used in making these figures are the corresponding values for each approximation (i.e. 
those shown in Table II) 



VI. CONCLUSION 



The theory of finite size scaling (FSS) has been demonstrated to be of great utility in analyzing numerical data 
of finite systems [1-8]. Moreover, it has been demonstrated that the theory also works for quantum systems, where 
the formalism applies to a complete basis expansion in the Hilbert space [24-33]. In this study we implemented the 
finite-element method which uses a local basis expansion consisting of arbitrary polynomials in real space and with 
variable resolution. It was demonstrated that this method can be used to compute matrix elements (O)^ (used 
in the finite size scaling formalism) accurately and with monotonic convergence with respect to the system size N. 
We have implemented the combined approach on different approximations to electronic calculations such as Hartree- 
Fock approximation (+ electron-electron correlation) and density-functional theory (DFT) under the local density 
approximation and found that even at the simplest level in DFT, the results we have obtained are in good agreement 
with other more accurate estimates such as large order perturbation theory [43]. 

The finite-element method has already shown great promise for extended systems since it uses local basis expansion 
to yield sparse matrices that lend the methodology easily to parallel architectures [14, 19]. Furtermore, the Khon-Sham 
equations in DFT provides a general method for extended system. Therefore we expect this combined approach can 
be applied to electron structure in atoms, molecules, and extended systems without loss of generality. Nevertheless, 
more studies need to be done to confirm the applicability of FSS with FEM for larger systems, where this combined 
approach can be used to shed some light into the nature of phase transitions. 
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